################################################################################################
#### 运行配置文件
################################################################################################

source /public/home/xxf2019/20231121_singleMuti/config/config.sh


################################################################################################
#### 结合RNA的细胞分类，计算atac按照501bp计算peak
################################################################################################

#### call peak
conda activate archr_v1
${Rscript_archr} ${scripts_dir}/run_archR.callPeak.R \
--rna_data_file ${input_dir}/testis_combined.Rdata \
--input_dir ${input_dir} \
--cluster_cell_file ${config_path}/cluster_celltype.csv \
--cpu 10 \
--out_path ${output_dir}/atac_res 

<<EOF
#### 基于atac的数据，构建每个细胞中peak的counts矩阵
for cluster in `cat ${config_path}/cluster_celltype.csv | awk -F, '{print $1}'`
do
echo " sh ${scripts_dir}/cluster_peak_count.sh ${cluster} ${config_path}" | \
qsub -N ${cluster}_"cluster_peak" -l nodes=1:ppn=30,mem=300gb,walltime=240:00:00 -q smp -d ${qsub_log_path}
done

#### 合并所有通过质控的细胞
cat ${output_dir}/cluster_all_result/cell_fragment/testis01_all_qc_barcode.csv \
> ${output_dir}/cluster_all_result/cell_fragment/testis_merge_all_qc_barcode_new.txt
cat ${output_dir}/cluster_all_result/cell_fragment/testis02_all_qc_barcode.csv | sed '1d' \
>> ${output_dir}/cluster_all_result/cell_fragment/testis_merge_all_qc_barcode_new.txt
cat ${output_dir}/cluster_all_result/cell_fragment/testis03_all_qc_barcode.csv | sed '1d' \
>> ${output_dir}/cluster_all_result/cell_fragment/testis_merge_all_qc_barcode_new.txt

#### 合并不同批次所有通过质控的细胞的counts的数据
for cluster in `cat ${config_path}/cluster_celltype.csv | awk -F, '{print $1}'`
do
echo " sh ${scripts_dir}/cluster_peak_countCombine.sh ${cluster} ${config_path}" | \
qsub -N ${cluster}_"countCombine" -l nodes=1:ppn=30,mem=300gb,walltime=240:00:00 -q student -d ${qsub_log_path}
done
EOF

################################################################################################
#### 标记并展示RNA聚类结果
################################################################################################
## 注释细胞类型
${Rscript_archr} ${scripts_dir}/annotation_celltype.R \
--input_file ${input_dir}/testis_combined.Rdata \
--cell_type_file ${output_dir}/atac_res/testis_merge_all_qc_barcode_new.txt \
--out_path ${input_dir} 

## 展示细胞亚群及marker基因
${Rscript_archr} ${scripts_dir}/rna_celltype_show.R \
--input_file ${input_dir}/testis_combined.annotationCellType.Rdata \
--out_path ${output_dir}/rna_celltype 

################################################################################################
## 质控ATAC细胞并重新call peak
## RNA的异常细胞也同时去掉
################################################################################################
## 去除异常的一簇细胞，该Patchytene只在testis01患者中有比例为11%（700多个细胞），远高于别的患者的2.5%（300多个），去除以后比例正常
${Rscript_archr} ${scripts_dir}/atac_cell_check.R \
--comine_data_file ${output_dir}/atac_res/testis_combined_peak.combineRNA.Rdata \
--rna_file ${input_dir}/testis_combined.annotationCellType.Rdata \
--out_path ${output_dir}/qc_atac
## ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata
## ${output_dir}/qc_atac/testis_combined.annotationCellType.qc.Rdata
## 20852个细胞

## 重新call peak和鉴定motif
${Rscript_archr} ${scripts_dir}/archr_callPeak.R \
--comine_data_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/qc_atac
## ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata

## 导出atac的矩阵
${Rscript_archr} ${scripts_dir}/get_genescorematrix.R \
--comine_data_all_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/qc_atac

## 展示细胞亚群及marker基因
## rna和atac均展示
${Rscript_archr} ${scripts_dir}/rna_celltype_show.R \
--rna_file ${output_dir}/qc_atac/testis_combined.annotationCellType.qc.Rdata \
--atac_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/celltype_plot


################################################################################################
#### 计算peak-gene的调控关系
################################################################################################
<<EOF
#### 填补细胞的atac矩阵
for cluster in `cat ${config_path}/cluster_celltype.csv | awk -F, '{print $1}'`
do
echo " sh ${scripts_dir}/cluster_peak_magic.sh ${cluster} ${config_path} " | \
qsub -N ${cluster}_"peak_magic" -l nodes=1:ppn=30,mem=300gb,walltime=240:00:00 -q student -d ${qsub_log_path}
done

#### 计算相关系数
for cluster in `cat ${config_path}/cluster_celltype.csv | awk -F, '{print $1}'`
do
echo " sh ${scripts_dir}/cluster_peak_correlation.sh ${cluster} ${config_path} " | \
qsub -N ${cluster}_"peak_correlation" -l nodes=1:ppn=30,mem=300gb,walltime=240:00:00 -q student -d ${qsub_log_path}
done



## 最后使用archR自带的
################################################################################################
#### 鉴定细胞特异的转录因子
################################################################################################

#### RG那一列替换为细胞类型所在的CB那一列
for bam_file in `ls ${input_dir} | grep sorted.bam | grep merge`
do
new_name=`echo ${bam_file} | awk -F"." '{print $1}'`
samtools view -h ${input_dir}/${bam_file} | \
awk 'BEGIN{OFS=FS="\t"} {for (i=1; i<=NF; i++) {if ($i ~ /^CB:Z:/) {gsub(/RG:Z:[^\t]+/, "RG:Z:"substr($i, 6)); break}} print}' | \
samtools view -bS - -o ${input_dir}/${new_name}.reviewRG.bam
samtools index ${input_dir}/${new_name}.reviewRG.bam
done

#### 提取每个细胞亚群特异的peak
mkdir -p ${output_dir}/motif
for cluster in ` cat ${config_path}/cluster_celltype.csv | awk -F',' '{print $1}' `
do
clustern=`echo ${cluster} | sed 's/cluster//'`
cat ${output_dir}/atac_res/testis_combined_peak.tsv | \
awk -F'\t' '{OFS="\t"}{ if($19==clustern) print $1,$2,$3 }' clustern=${clustern} > ${output_dir}/motif/${cluster}_peak.bed
done

#### 基于单细胞的ATAC数据，需要提供peak和bam，每个细胞类型分别鉴定peak对应的motif以及motif在每个细胞的活性程度z-score
## 参考以下内容写代码
## https://greenleaflab.github.io/chromVAR/articles/Articles/Applications.html#differential-accessibility-and-variability
## 文章：https://www.nature.com/articles/nmeth.4401
genome_type="hg38"
cpu=10

for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' | grep -E "cluster9|cluster15|cluster3|cluster2" `
do
clus=`echo ${line} | awk -F, '{print $1}'`
celltype_name=`echo ${line} | awk -F, '{print $2}'`

echo ${line}
echo " sh ${scripts_dir}/run_chromvar.sh ${clus} ${celltype_name} ${genome_type} ${cpu} ${config_path} " | \
qsub -N ${clus}_"chromvar" -l nodes=1:ppn=20,mem=200gb,walltime=240:00:00 -q smp -d ${qsub_log_path}
done

#### 1、计算TF的活性和基因表达的相关性
#### 2、根据peak和基因的关系以及peak和TF的关系，计算TF在该基因的激活程度，TF的z-score
cor_gene_peak=0.5
cpu=5

for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
	
cluster=`echo ${line} | awk -F, '{print $1}'`
celltype_name=`echo ${line} | awk -F, '{print $2}'`

echo ${line}
echo " sh ${scripts_dir}/tf_gene.sh ${cluster} ${cor_gene_peak} ${celltype_name} ${cpu} ${config_path} " | \
qsub -N ${cluster}_"tf_gene" -l nodes=1:ppn=16,mem=200gb,walltime=240:00:00 -q smp -d ${qsub_log_path}

done
EOF



################################################################################################
#### 拟时序
################################################################################################
<<EOF
## 基于表达数据建立拟时序
${Rscript_archr} ${scripts_dir}/run_monocle.R \
--input_file ${input_dir}/testis_combined.annotationCellType.Rdata \
--out_path ${output_dir}/monocole

## 整合基因的atac、rna以及拟时序的结果为Rdata
${Rscript_archr} ${scripts_dir}/run_combine_rna-atac.addPseudotime.R \
--rna_data_file ${input_dir}/testis_combined.annotationCellType.Rdata \
--rna_data_monocle_file ${output_dir}/monocole/testis.monocle.Rdata \
--atac_data_file ${output_dir}/atac_res/testis_combined_peak.Rdata \
--time_diff_gene_file ${output_dir}/monocole/pseudotime_ordergene.tsv \
--gff3_file ~/ref/GTF/gencode.v32.annotation.gff3 \
--gene_ensg_file ~/ref/GTF/gencode.v32.gene_ensg.tsv \
--out_path ${output_dir}/atac_res

## rna拟时序画图
## RNA本身的以及，与时序相关基因的表达和峰开放
${Rscript_archr} ${scripts_dir}/run_monocle_plot.R \
--rna_data_monocle_file ${output_dir}/monocole/testis.monocle.Rdata \
--comine_data_file ${output_dir}/atac_res/testis_combined_peak.combineRNA.Rdata \
--time_diff_gene_file ${output_dir}/monocole/pseudotime_ordergene.tsv \
--out_path ${output_dir}/monocole
EOF

## 基于ATAC的时序分析
## 最后使用
${Rscript_archr} ${scripts_dir}/run_trajectory.R \
--comine_data_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata \
--out_path ${output_dir}/trajectory


################################################################################################
## 分每个细胞亚群,分别鉴定motif以及peak2gene以及TF-gene的对应关系
################################################################################################
##########################################
## 鉴定motif以及peak2gene
## peak2gene默认是250kb内
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
echo " sh ${scripts_dir}/chromvar_p2gene.sh ${cluster} ${config_path} " | \
qsub -N ${cluster}_"chromvar_p2gene" -l nodes=1:ppn=10,mem=100gb,walltime=240:00:00 -q batch -d ${qsub_log_path}
done

##########################################
## 鉴定TF-gene的可能关系
## 默认代码
## TF选择在所有细胞类型间存在差异的
## 先跑,结果比较结果后面会进行调整
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
echo " sh ${scripts_dir}/run_TFregulators.sh ${cluster} ${config_path} " | \
qsub -N ${cluster}_"run_TFregulators" -l nodes=1:ppn=10,mem=100gb,walltime=240:00:00 -q batch -d ${qsub_log_path}
done

<<EOF
##########################################
## 画所有的TF-gene，对TF不筛选
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
echo " sh ${scripts_dir}/run_TFregulators.All.sh ${cluster} ${config_path} " | \
qsub -N ${cluster}_"run_TFregulators" -l nodes=1:ppn=10,mem=100gb,walltime=240:00:00 -q batch -d ${qsub_log_path}
done
EOF

##########################################
## 用新的参数鉴定TF-gene的可能关系
## 更宽松，能鉴定更多已报道的可靠的TF
cor=0.3
maxDelta=0.7
for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
echo " sh ${scripts_dir}/run_TFregulators.use.sh ${cluster} ${config_path} ${cor} ${maxDelta}" | \
qsub -N ${cluster}_"run_TFregulators" -l nodes=1:ppn=10,mem=100gb,walltime=240:00:00 -q smp -d ${qsub_log_path}
done

## 所有细胞类型中阳性的TF合并到一张表格
## 头行
cat ${output_dir}/tf_regulators_${cor}_${maxDelta}/cluster0/TFregulatorPlots_All/*tsv | head -1 | \
awk -F'\t' '{OFS="\t"}{print "TF",$0,"cluster","cell_type"}' \
> ${output_dir}/tf_regulators_${cor}_${maxDelta}/Positve_TF-Gene.tsv

for line in ` cat ${config_path}/cluster_celltype.csv | tr ' ' '_' | tr '&' '-' `
do
cluster=`echo ${line} | awk -F, '{print $1}'`
cell_type=`echo ${line} | awk -F, '{print $2}'`
echo ${cell_type}
for file in `ls ${output_dir}/tf_regulators_${cor}_${maxDelta}/${cluster}/TFregulatorPlots_All/ | grep tsv`
do
tf=`echo ${file} | awk -F'_' '{print $1}'`
cat ${output_dir}/tf_regulators_${cor}_${maxDelta}/${cluster}/TFregulatorPlots_All/${file} | grep -v putative_targets | \
awk -F'\t' '{OFS="\t"}{print TF,$0,cluster,cell_type}' TF=${tf} cluster=${cluster} cell_type=${cell_type} \
>> ${output_dir}/tf_regulators_${cor}_${maxDelta}/Positve_TF-Gene.tsv
done
done
## 提取存在调控的
cat ${output_dir}/tf_regulators_${cor}_${maxDelta}/Positve_TF-Gene.tsv | grep -E "YES|putative_targets" \
> ${output_dir}/tf_regulators_${cor}_${maxDelta}/Positve_TF-Gene.onlyTargetGene.tsv

##########################################
## 所有已报道的关键的TF在motif水平、其本身开放程度以及表达水平进行染色
${Rscript_archr} ${scripts_dir}/plotEmbedding_TFregulators.R \
--rna_file ${input_dir}/testis_combined.annotationCellType.Rdata \
--comine_data_all_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata \
--report_tf_file ${config_path}/TF_upregulated_during_spermatogenesis.csv \
--out_path ${output_dir}/report_tf

################################################################################################
#### 标记并展示已知的重要基因在rna及atac的热图
################################################################################################
## 基因分为总的，TF以及nonTF
## 展示按照基因的表达聚类以及按照基因的固定顺序
${Rscript_archr} ${scripts_dir}/pheatmap_knowngene.R \
--rna_file ${input_dir}/testis_combined.annotationCellType.Rdata \
--comine_data_all_file ${output_dir}/qc_atac/testis_combined_peak.combineRNA.qc.Rdata \
--known_gene_file ${config_path}/SCOS_known_genes.csv \
--out_path ${output_dir}/pheatmap_knowngene


################################################################################################
## 共享数据
mkdir -p ${output_dir}/shareData/peak-gene
mkdir -p ${output_dir}/shareData/motif/cisbp
mkdir -p ${output_dir}/shareData/motif/jaspar

cd ${output_dir}/shareData

## 细胞类型注释
cp -rf ${config_path}/cluster_celltype.csv ${output_dir}/shareData/

## 转录组数据
ln -snf ${input_dir}/testis_combined.annotationCellType.Rdata ${output_dir}/shareData/

## atac包含peak的信息
ln -snf ${output_dir}/atac_res/testis_combined_peak.Rdata ${output_dir}/shareData/

## peak-gene
for file in `find ${output_dir}/cluster_all_result | grep correlation_magic.Rdata`
do
ln -snf ${file} ${output_dir}/shareData/peak-gene
done

## motif
for file in `find ${output_dir}/motif | grep rda | grep cisbp`
do
ln -snf ${file} ${output_dir}/shareData/motif/cisbp
done
for file in `find ${output_dir}/motif | grep rda | grep jaspar`
do
ln -snf ${file} ${output_dir}/shareData/motif/jaspar
done